The Hubbard model with intersite interaction 
within the Composite Operator Method 
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We study the one- and two- dimensional extended Hubbard model by means of the Composite 
Operator Method within the 2-pole approximation. The fermionic propagator is computed fully 
self-consistently as a function of temperature, filling and Coulomb interactions. The behaviors of 
the chemical potential (global indicator) and of the double occupancy and nearest-neighbor density- 
density correlator (local indicators) are analyzed in detail as primary sources of information regard- 
ing the instability of the paramagnetic (metal and insulator) phase towards charge ordering driven 
by the intersite Coulomb interaction. Very rich phase diagrams (multiple first and second order 
phase transitions, critical points, reentrant behavior) have been found and discussed with respect to 
both metal-insulator and charge ordering transitions: the connections with the experimental find- 
ings relative to some manganese compounds are analyzed. Moreover, the possibility of improving 
the capability of describing cuprates with respect to the simple Hubbard model is discussed through 
the analysis of the Fermi surface and density of states features. We also report about the specific 
heat behavior in presence of the intersite interaction and the appearance of crossing points. 



I. INTRODUCTION 

Many authors have emphasized the importance of 
considering non-local Coulomb interactions in describ- 
ing doped systems like cuprate superconductors or 
fuUeridesiiSiSi^. The simplest Hamiltonian satisfying 
these requirements is the extended Hubbard model, 
where a nearest-neighbor Coulomb interaction term is 
added to the original Hubbard Hamiltonian®. The in- 
clusion of non-local Coulomb interactions substantially 
modifies the electronic properties of the model. For in- 
stance, the charge transfer excitons, which can only be 
detected by optical spectroscopy at half filling, attain 
some charge in doped systems and become visible in di- 
rect and inverse photoelectron spectroscopies^. Other 
studies, using an effective extended Hubbard model, sup- 
port the appearance, upon doping, of states evenly dis- 
tributed inside the gap*. This suggests that the general 
features of the cuprates can be well described by using 
this effective Hubbard Hamiltonian, which has also been 
used to mimic some of their experimental features in the 
superconducting state by means of a BCS treatmenljS. 

The Hubbard model with intersite Coulomb interac- 
tion is also one of the simplest models capable to de- 
scribe charge ordering (CO) in interacting electron sys- 
tems. Already in 1938 Wigner proposediS that a low 
density interacting electron gas crystallizes in a lattice 
in order to minimize the Coulomb repulsion. At higher 
densities crystallization is possible if the kinetics energy 
is reduced by spin or phonon interactionsii. Charge 
ordering has been experimentally observed in a vari- 
ety of systems: GaAs/AlGaAs heterostructuresi^, rare- 
earth pnictides like Yb4^Asj^, colossal magnetoresistance 
compoundsi^, unconventional spin-Peierls materials a — 
NaV205^^ , cuprates^®, manganitesii, magnetitei^, vana- 
dium oxidesiS, Bechgaard salts^S. 

Among the many analytical methods used to study 



the extended Hubbard model we recall: Hartree- 
Fock approximations!, perturbation theorjiSS, dynamical 
mean field theorjs^S,, slave boson approach2^*2^, coherent 
potential approximation^®. Numerical studies by means 
of Quantum Monte CarloS, Lanczos technique^^ and ex- 
act diagonalizationS^ have also to be recalled. 

In this manuscript, we analyze the extended Hubbard 
model by means of the Composite Operator Method 
(COM) within the 2-pole approximation (see Refs. lSClSlI 
and references therein). 

The COM is based on two main concepts: (i) the exci- 
tations present in interacting systems are far from being 
the original electrons and need to be described by asymp- 
totic fields in the form of composite operators; (ii) the use 
of composite operators requires the enforcement of the 
non-canonical algebra they obey in order to properly fix 
the representation where their propagators are realized. 
This latter task is effectively undertaken by computing 
the values of the unknown correlators appearing in the 
calculations by means of Algebra Constraints^. 

The detailed analysis of the instabilities of the ho- 
mogenous paramagnetic phase of the extended Hubbard 
model towards charge ordered inhomogeneous phases is 
the main task undertaken with this manuscript. With 
respect to this, we have analyzed the rank of the tran- 
sitions and their relations with the metal-insulator one. 
Actually, this study is the only relevant at room tem- 
peratures and/or in presence of frustration as any spin 
ordered phase (e.g., the antiferromagnetic phase) would 
be inaccessible in such situations. It is worth noticing 
that the very rich phase diagrams (multiple first and 
second order phase transitions, critical points, reentrant 
behavior) that have been found can be put in connec- 
tion with the experimental findings relative to some man- 
ganese compounds'^'^'^'^''^^. We have also discussed in de- 
tail the possibility of improving the capability of describ- 
ing cuprates with respect to the simple Hubbard model 
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and the appearance of crossing points in the specific heat 
in presence of the intersite interaction. 

In the next section we present the model, the basis 
and the solution according to the COM. In the subse- 
quent sections, we comment our results for the chemical 
potential, the phase diagram, the double occupancy, the 
nearest-neighbor density-density correlator, the internal 
and kinetic energies, the Fermi surface, the density of 
states and the specific heat. The types of charge or- 
dered phases found according to the sign of the intersite 
potential are described, the rank of the transitions is ev- 
idenced through the analysis of the discontinuities in the 
chemical potential, double occupancy, nearest-neighbor 
density-density correlator and kinetic and internal ener- 
gies, quite complex phase diagram (with metal to insu- 
lator and to charge ordered phase transitions and reen- 
trant behavior) are drawn and commented, the behavior 
of single- (double occupancy) and two- (nearest-neighbor 
charge) site correlators is studied in order to get infor- 
mation about the actual charge distribution, the value of 
the filling at which the nesting appear is determined as a 
function of intersite potential, the specific heat behavior 
is studied in comparison to that of the simple Hubbard 
model. 



of motion of this latter take the general form 



(2.3) 



where e(k) = ^[e(i,j)] is the n x n energy matrix de- 
scribing the projected dynamics. The energy matrix 
can be computed as e(k) = TO(k)/^^(k) where /(k) = 
J-{{ip{i,t),ip^{i,t)}) is the normalization matrix of the 
basis and m(k) = T{{i^^{i,t),^\i,t)}). If we neglect 
Sj{i) we obtain a pole structure for the retarded thermal 
Green's function G(k,tj) = JP(7^[?^(^)?^^(j)]) {T^ is the 
retarded operator) that will obey the following equation 
of motion 



ujG{k, Lu) = /(k) + e(k)G'(k, u) 
The solution of Eq. (|2.4|l is 

(7W(k) 



^ ^ Z,? — 



K.(k) 



(2.4) 



(2.5) 



where Ei(k) are the eigenvalues of e(k) and the spectral 
weights cr^'^ (k) can be computed as 



II. HAMILTONIAN, FIELD EQUATIONS AND 
SOLUTION 



We will study a generalized version of the Hub- 
bard model^ which includes the intersite Coulomb 
interactio n?^'?i??-?7i??i??i49-4i . Accordingly, the Hamilto- 
nian under analysis reads as 



H = ^[-^ct(i)c(i) - 2dtc\i)c°'{i)] 
+ ^[f/nt(i)ni(i)+dFn(i)n"(z)] 



(2.1) 



where /i is the chemical potential, cJ{i) = (c|(i)c|(i)) is 
the creation electronic operator in spinorial notation, i — 
(i,t), i is one lattice vector of the d-dimensional square 
lattice, t, as usually done in the related literature, is both 
the time variable and the hopping integral, the context 
will clarify the use, U is the onsite Coulomb interaction, 
V is the intersite interaction, n{i) = n|(«) -I- ?t-j,(«)j 't-ct(*) 
is the number operator for electrons of spin a. Hereafter, 
t will be used as reference unit for all energies. We have 
used the notation 



6"(i,t)-^aij<^(j,i) 



(2.2) 



where (j) can be any operator and ay is the projector 
on the first 2d neighbor sites on the lattice. We have 
a(k) = .?"[aij] = cos(fc„), where T is the 

Fourier transform. 

Within the Composite Operator Method2&, once we 
choose a n-component spinorial basis '0(*): the equations 



'-i2(k) 



E 

c=l 



0„(k)f)-i(k)/,b(k) a, 6=1, 



(2.6) 



where the matrix ri(k) has the eigenvectors of e(k) as 
columns^. In this manuscript, we will use a 2-pole ap- 
proximation within the COM. The reader interested to 
more elaborate self-energy treatments (in order to take 
into account 5j{i) we should compute the higher order 
propag ator ('R\5 i(i)5i'' (i)])) should refer, for instance, 
to R,ef. l42l48l44l 

For the model under analysis in this manuscript, we 
choose, as basic fields, the Hubbard operators, which 
faithfully describe the Hubbard subbands as eigenopera- 
tors of the ionic model, 



7j{i) 



(2.7) 



where — n{i)c{i) and r/{i) = c{i) — ^(i) = [1 — 
n{i)]c{i). They satisfy the following equations of motion 



-fi^i) - 2d[tc«(z) + t7r{i) - Vai)n"{i)]\ 
-{n-U)T]{i)+2d[tTT{i) + Vr]{i)n°'{i)] I ^^''^^ 



where 



n{i)^laf'n^{i)c"{i)+c{i)c''\i)c{i) 



(2.9) 



a'^ = (— 1,(t), np(i) = {n{i),n{i)) is the charge and spin 
number operator, n{i) — c^{i)ac{i) and a are the Pauli 
matrices. 
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After the choice we made for the basis, we can com- 
pute, in the 2-pole approximation within COM, the en- 
ergy spectra Ei(k) and the spectral density functions 
CT'^*-'(k) according to the general procedure given above. 
The lengthy expressions can be found in Appendix. It 
is worth noting that some parameters, not connected to 
the Green's function under analysis, appear in the ex- 
pressions: Xc — {n{i)n°'{i)) and p = j{n"{i)nf^{i)) — 

([c|(i)c|(i)]"c|(i)c|(i)). The first one, x": "^ill be com- 
puted by calculating the density-density correlation func- 
tion {n(i)n{j)) within the one-loop approximation"*^ (see 
second equation in Eqs. (jAlOp ). The second will be fixed 
by the local algebra constraintSS (^(i)?7^(i)) = 0. Such 
a procedure is very peculiar to the COM treatmen1»22. as 
discussed in the introduction. By solving the set of cou- 
pled self-consistent equations (|A9IA10|I , we can calculate 
the various correlation functions and the physical prop- 
erties of the system. Results will be presented in the 
following sections for the one-dimensional (ID) and the 
two-dimensional (2D) systems. 

Before moving to the results, it is worth noting that 
the set of self-consistent equations (jA9IA10|) is highly 
non-linear. According to this, it is natural to expect a 
certain number of coexisting solutions with quite differ- 
ent features. Actually, the set admits only two distinct 
solutions that, hereafter, we will call COMl and COM2, 
according to the main sign of parameter p. In partic- 
ular, as a function of the filling, we have a p positive 
and of the order of the filling in COMl and a p negative 
or very small and positive in COM2. As the Compos- 
ite Operator Method tries to give answers in the whole 
space of model and physical parameters and as the Hub- 
bard model response is profoundly different according to 
the region of the latter space, the presence of two so- 
lutions, so different in their features, should be seen as 
a richness of the method. Due to the difficulties inher- 
ent to the task of studying the whole phase diagram, 
many other approximations focus just on one region and 
usually give wrong results in the rest of the parameter 
space. COM, also thanks to the opportunity of choos- 
ing between the two solutions according to the features 
one expects to describe, has proved to be capable to ex- 
plore the whole phase diagram, to get significatively good 
results for all dimensions and values of the tuning param- 
eters, to give relevant interpretations of the experimental 
facts and to constitute a solid basis for further investi- 
gations. It is also worth remembering that, at the end 
of the day, our analysis has, as main goal, the descrip- 
tion of features, more or less anomalous, of real physical 
systems by means of Hamiltonian models. We can never 
be sure that our Hamiltonian contains all, and only, rel- 
evant ingredients. Different experimental situations can 
be described by different analytical (and/or numerical) 
solutions of the same Hamiltonian. For a fixed set of 
parameters (filling, temperature, pressure, ...), the dif- 
ferent solutions will have different free energies and one 
can simply think to choose one or the other according to 
this. This is the correct procedure if we wish to answer 



the question: "Which is the phase effectively described by 
this Hamiltonian under such external conditions?" . On 
the other hand, a solution different from the one with 
lower free energy can better describe the real experimen- 
tal situation. This latter, although mainly determined 
by the ingredients already present in the Hamiltonian 
(the possibility to find a solution with similar character- 
istics assures this), is that effectively realized in nature 
owing to some marginal interactions not present in the 
Hamiltonian chosen. If we would include such interac- 
tions in the Hamiltonian and compute new solutions, we 
will find that the lower free energy one will be the one 
capable to describe the actual physical situation. Now, 
the most important consideration to be done is that this 
latter solution will be essential identical to that obtained 
with the previous Hamiltonian and that was already ca- 
pable to describe the actual physical situation although 
it has not the lower free energy. According to this, in 
the comparison with the experiments no solution can be 
discarded a priori (neither on the basis of the free energy 
determination) and, with this idea in mind, we have here, 
and in many other works, presented the results for the 
two solutions we have obtained. 

For the one-dimensional system we will consider only 
COM2 solution as our past experience suggests that this 
is the one best suited to describe the physics of the 
simple one-dimensional Hubbard model for which we 
got excellent agreements with the Bethe Ansatz exact 
solution^SiiSi^SJ^ . For the two-dimensional case, we will 
study both solutions as they will permit us to analyze, in 
particular in proximity of half filling and as regards the 
metal-insulator transition and its relation to the charge 
ordering transition, two different behaviors that could 
be both observed experimentally. Also in this case the 
many positive comparisons with numerical results, that 
we have obtained in the previous studies on the sim- 
ple two-dimensional Hubbard model'^^i^°i^^i^^i^^i^'*, will 
be used as a guide throughout all the analysis. 



III. THE CHEMICAL POTENTIAL 

The chemical potential can be determined as a function 
of the parameters n, T, U , V by solving the system of 
self-consistent equations (|A9IA10|I . It is worth noticing 
that our results show that the relation 

H{2 - n) ^ U + AdV - fi{n) (3.1) 

required by the particle-hole symmetry, is exactly sat- 
isfied. In particular at half filling we have fi{l) — 
+ 2dV. This is due to the fact that among the possi- 
ble representations the Algebra Constraint coming from 
the Pauli principle [see first equation in Eqs. (|A10|) ] se- 
lects the one which preserves the particle-hole symmetry 
of the model-°. Any other choice for the equation fixing 
the parameter p leads to a violation of the symmetry. 
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FIG. 1: The chemical potential at [/ = 2 as a function of: 
(top) the particle density for V — 0, 1; (bottom) the intersite 
Coulomb interaction for several values of n. 



A. One-dimensional system 

The chemical potential as a function of the particle 
density is reported in Fig. ^ (top) for J7 = 2 and V = 
and 1. For any value of the filling the chemical potential 
increases by increasing the intersite Coulomb interaction 
and the increment is an increasing function of the filling: 
within the paramagnetic phase the average probability 
for two particles to be nearest neighbors increases with 
the filling and accordingly increases the free energy in 
the presence of a repulsive intersite Coulomb interaction, 
then, the behavior of the chemical potential follows. 

For negative (i.e., attractive) values of the intersite 
Coulomb potential, we can see that the paramagnetic 
solution is unstable towards a phase with charge separa- 
tion (i.e., with charge ordering). As a matter of fact, an 
attractive potential between charges at nearest neighbor 
sites favors a rearrangement of the particles in a one- 
particle-per-site scheme in order to maximize the gain in 
energy. At fillings less than one this scheme can lower the 
energy more and more by accepting more particle in the 
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FIG. 2: The discontinuity of the chemical potential at half 
filling and T — 0.01 is plotted versus the intersite Coulomb 
interaction for U — 2, 3 and 4. 



system and increasing the number of occupied couples of 
nearest neighbor sites. According to this, the chemical 
potential shows a negative slope which, on increasing the 
value of the intersite potential, manifests at lower and 
lower values of the filling. The critical value of the inter- 
site Coulomb potential which controls the transition to 
the charge ordered state depends on the intensity of the 
local Coulomb interaction. In order to illustrate the case, 
in Fig. n (bottom), the chemical potential is given versus 
intersite Coulomb potential for [/ = 2 and various values 
of the particle density: all curves cross at a certain value 
V « —0.5, below which the system exhibits a negative 
compressibility. 

As expected from the exact Bethe Ansatz solution of 
the one-dimensional simple Hubbard model, for any U > 
there is a discontinuity at half filling, signalling a gap in 
the density of states and the occurrence of an insulating 
phase. As it can be seen from Fig. ^ (top), the effect 
of the intersite term V is to reduce the size of the gap. 
This is studied in Fig. |21 where the discontinuity of the 
chemical potential at half filling 



A/i = M+(l)-M-(l) 



(3.2) 



is plotted versus the intersite Coulomb potential for var- 
ious values of the onsite Coulomb potential at half filling 
and T = 0.01. We see that for a given value of the 
onsite Coulomb potential, there is a critical value Vc- 
For values of intersite Coulomb potential greater than 
Vc, the paramagnetic insulating phase becomes unstable 
(the chemical potential gets a negative slope) and there 
is a phase transition to a charge ordered insulating state 
(CO)^^. This kind of ordering is much different than 
that discussed previously (i.e., the one-particle-per-site 
type); in this case the repulsion among particles favors 
a checkerboard pattern with half sites double occupied 
(see in the next sections the discussion about the double 
occupancy) and half empty, these latter ones being the 
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FIG. 3: The discontinuity of the chemical potential at half 
filling is plotted versus the intersite Coulomb interaction for 
various values of U and T = 0.01. 



nearest neighbors of the former ones. By analyzing the 
derivative [d{AiJ,)/dV]v=Va^ we can show that at zero 
temperature the transition is second order for U < 2 and 
first order for higher values of the onsite Coulomb poten- 
tial. 



B. Two-dimensional system 




As already noted for the ID system, also in 2D within 
COM2 solution the discontinuity A/x of the chemical po- 
tential at half filling decreases by increasing the value of 
the intersite Coulomb potential and the system exhibits 
a phase transition to a charge ordered phase at some crit- 
ical value Vc- The values of Vc, as a function of the onsite 
Coulomb potential, are shown in Fig. 01 where A/i is re- 
ported versus the intersite Coulomb potential for various 
values of the onsite interaction U. 

Within COMl solution, the behavior of the chemical 
potential as a function of the particle density is shown in 
Fig. ^ (top) for [/ = 4 and various values of the intersite 
Coulomb potential. We see that for attractive values of 
the intersite Coulomb interaction the chemical potential 
decreases by increasing the filling, showing an instabil- 
ity of the paramagnetic case towards phase separation. 
In particular, exactly at quarter filling (n=0.5), we can 
have a charge-ordered state of one-particle-per-site type. 
Away from quarter filling, the phase separation is be- 
tween two phases with different particle densities (these 
latter should be determined through a Maxwell construc- 
tion), whose nature could be investigated (they can be 
both charge-disordered or one ordered and another disor- 
dered) within a treatment where translational invariance 
should be relaxed. The instability of the paramagnetic 
phase can be here studied by plotting the chemical po- 
tential as a function of the intersite Coulomb potential. 
This is shown in Fig. 0] (bottom) where we see that the 



FIG. 4: The chemical potential at T = 0.01 and t/ = 4 versus: 
(top) the particle density for several values of V; (bottom) the 
intersite interaction for several values of n. 
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FIG. 5: The chemical potential versus the particle density at 
T = 0.01, V = 0.5 and several values of U. 
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curves for different values of filling all cross at some criti- 
cal value V* , depending on the onsite Coulomb potential 
[V* w -1.1 for U = 4]. For \V\ < \V*\ the chemical 
potential is a decreasing function of the filling. 

As in the simple Hubbard niodel45-50.5i,52,53,54^ COMl 
solution describes a metallic phase in the low regime of 
on-site Coulomb repulsion and a transition to the insulat- 
ing state when the potential reaches some critical value. 
This feature survives also when the intersite Coulomb in- 
teraction is taken into account. For a fixed value of this 
latter, there is a critical value of the on-site repulsion 
such that the chemical potential exhibits a discontinuity 
at half filling. This latter signals the opening of a gap 
in the density of states and therefore a Mott transition. 
This feature is illustrated in Fig. |31 Now, the current 
study shows that the critical value Uc decreases by in- 
creasing the intersite potential V. By further increasing 
the intersite potential, the system undergoes a second 
transition to a charge ordered state of checkerboard type. 
This transition is characterized by a discontinuity in the 
double occupancy as commented above and shown in the 
next sections. 
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FIG. 6: The critical value Vc where there is a phase transition 
from the Mott insulator to charge ordered state is plotted as 
a function of the onsite Coulomb potential at half filling and 
zero temperature. 



IV. PHASE DIAGRAMS 

A. One-dimensional system 

As reported in the previous section, COM2 solution 
for the one dimensional system shows that there is phase 
transition from the Mott insulating phase to a inhomo- 
geneous charge ordered state of checkerboard type for 
positive values of the intersite Coulomb potential greater 
than some V^- This result is consistent with many other 
studieo^^i^'^i^^^^i^^i'^^i^^i^". The nature of this phase tran- 
sition is not well understood yet and currently under in- 
tense investigatiorSiSiSiSSiSi. The phase diagram in the 
plane V-U is shown in Fig.|Hl where the critical value Vc 
of the intersite Coulomb potential is plotted as a function 
of the onsite potential [/. The arrow indicates the point 
where the phase transition from second order changes to 
first order. 

It is necessary noticing that this transition is quite 
different from that usually observed in proximity of the 
[/ = 21^ line between an inhomogeneous charge ordered 
phase and an homogenous spin ordered phase. In this 
manuscript, we decided to focus on the homogenous para- 
magnetic phase, its instabilities towards charge ordered 
inhomogeneous phases (paying attention to the rank of 
the transition) and the relation between these latter and 
the metal-insulator transition. These kinds of transitions 
are the only ones relevant at room temperatures and/or 
in presence of frustration as the antiferromagnetic phase 
is quite depressed in such cases. 

In Fig. Ul we give the phase diagram in the plane T- 
V ioT U — 2 and [/ = 8. For U = 2 the transition is 
second order for all values of T. For [/ = 8 the tran- 
sition is first order for T < 0.55 and second order for 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

1.0 
0.8 
0.6 



Mott insulator 








CO 








ID 






COM2 






n=l 


~ it 




U=2 



2.4 2.6 



0.0 



10 



2.8 3 
V 



3.2 3.4 



Mott insulator 






CO - 








ID 




COM2 




n=l 




U=8 



10.1 



10.2 10.3 
V 



10.4 



FIG. 7: The critical temperature Tc for the Mott insulator to 
charge ordered state phase transition is plotted as a function 
of the intersite Coulomb potential at half filling and U — 2 
(top) and U = 8 (bottom). 
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FIG. 8: The nearest-neighbor density-density correlation 
function x" is plotted as a function of V/Vc at half filling, 
T = 0.01 and (7 = 2, 8 and 20. 

T > 0.6. The arrow indicates the temperature where, 
by increasing T, the transition becomes continuous. It 
is interesting to observe that the transition is continuous 
if there is no reentrant temperature behavior. By this 
latter we mean a situation in which by increasing tem- 
perature we can first enter and then exit a phase when 
within another (e.g., see Fig. (bottom) at V = 10.1 
for increasing temperature from zero). When there is a 
reentrant temperature behavior, the transition is discon- 
tinuous up to the turning point, then becomes continu- 
ous. A reentrant temperature behavior has been experi- 
mentally observed^SiSiS^ and will be further discussed for 
the 2D system. It is worth noting that the transition is 
also marked by a discontinuity in the nearest-neighbor 
density-density correlation function Xc ~ {n" {i)n(i)) . 
This quantity is shown in Fig. |S1 as a function of V/Vc 
at half fiUing, T = 0.01 and C/ = 2, 8 and 20. At the 
transition Xc is continuous for U = 2 and discontinuous 
for U > 2. The nearest- neighbor density-density corre- 
lation function is a decreasing function of the intersite 
potential as the repulsion diminishes the probability of 
finding neighboring sites occupied. In particular, at the 
transition this probability reduces with a change of con- 
cavity (second order transition) or with a discontinuity 
(first order transition). 



B. Two-dimensional system: COM2 

The COM2 solution for the two-dimensional case has 
similar characteristic to the ID case. The phase diagrams 
in the planes V-U and T-V are shown in Figs. El (top) 
and 1^1 (bottom), respectively. 

At zero temperature the transition is continuous for 
U < 1.8 and first order for U > 1.9. For finite tem- 
perature and U — 8 a, reentrant behavior as func- 
tion of temperature is observed. The transition is 
first order up to the turning point T = 0.6, then be- 
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FIG. 9: (top) The critical value Vc for the Mott insulator to 
charge ordered state phase transition is plotted for the two- 
dimensional case (COM2 solution) as a function of: (top) the 
onsite Coulomb potential at half filling and zero temperature 
(COM2 solution); (bottom) V at half filling and U = 8. 



comes continuous. For U < 1.8 no re-entrant behav- 
ior is observed. The fact that charge ordering may 
disappear by decreasing temperature has been experi- 
mentally observed in Pro. 65(Cao.7S'ro. 3)0.35^^03^ and 
La2-2xSri+2xMn20r (0.47 < a; < 0.62)^^. On the 
theoretical side a re-entrant temperature has been ob- 
tained in the context of the extended Hubbard model at 
quarter filling^'^-^^-^'*. 

In Fig. ^1 (top) the nearest-neighbor density-density 
correlation function Xc — ('^"(*)'^(*)) plotted as a func- 
tion of V/Vc at half filling, U = 8 and various tempera- 
tures. At the transition Xc discontinuous for T ~ 0.01 
and T = 0.4, and continuous for T — 1, in agreement 
with the phase diagram shown in Fig. |51 (bottom). 



C. Two-dimensional system: COMl 

The results given for the chemical potential show that 
for low values of the on-site Coulomb interaction the sys- 
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FIG. 10: (top) The nearest-neighbor density-density correla- 
tion function is plotted as a function of V jVc at half filling, 
f/ = 8 and T = 0.01, 0.4 and T = 1, for the two-dimensional 
case (COM2 solution), (bottom) The energy gap Ai? at half 
filling is reported as a function of U for several values of U for 
T=0, for the two-dimensional case (COMl solution). 



tern is in a metallic state and undergoes a metal-insulator 
transition for a critical value C/c, which depends on the 
intensity of the intersite potential. In order to study this 
metal-insulator transition (MIT) we consider the quan- 
tity 



l\E ^ £'i(0,0) - £'2(7r,7r) 



(4.1) 



which measures the distance between the bottom of the 
upper subbands and the top of the lower one. After 
Eq. HA4p . at half filling we have 



i?(k) = i?ia(k) 

= go 



mi2(k) = TOiQ:(k) m 



Therefore, we simply have 



i?i ^ -Stp + W{C^^ + Ci"i) 

.go = -t/ + 8ni-X?) (4.2) 
1 



AE = 2Ri 
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FIG. 11: The complete phase diagram in the plane U-V 
is shown at half filling for zero temperature, for the two- 
dimensional case (COMl solution). 



For a fixed value of the onsite Coulomb repulsion, 
the equation AE = will determine the critical value 
Vc where a metal-insulator transition occurs. The re- 
sults show that the metallic region is compressed by the 
presence of the intersite interaction and disappears for 
V > 5.7. This is shown in Fig. ^1 (bottom), where AE is 
reported versus the onsite Coulomb potential for various 
values of the intersite Coulomb potential at zero temper- 
ature. 

For V > h.7 the equation AE ~ does not have a 
solution and the metallic phase disappears. By further 
increasing the intersite Coulomb potential, there is an 
instability in the self-consistent equations ljA9IA10|) . the 
double occupancy exhibits a discontinuity leading the 
system to a charge ordered state. The complete phase 
diagram in the plane V-U is shown in Fig. ^2 The di- 
agram is characterized by two critical curves, Vd and 
Vc2, which separate the different phases. Vd controls the 
MIT and Vc2 controls the transition to a charge ordered 
state. At Vci the transition is first order for U < 12 and 
second order for U > 12.2; at Vc2 we have first order 
for U > 1.9 and second order for U < 1.8. As proposed 
in RefT lG^ the insulating state can be characterized by 
the order parameters {^aij)^Ujodd)) and {r]aU)vlijodd)) 
which vanish at the MIT (jodd is any site reachable in an 
odd number of hops from site j). The nearest neighbor 
hopping ampHtudes A — = iv" and 

B = {ri°'{i)^Hi)) = {£."{i)vH^)) are reported in Figs. [El 
as functions of the onsite Coulomb potential for various 
values of the intersite Coulomb potential. We see that 
for a given value of the intersite Coulomb potential, the 
probability amplitude A suddenly vanishes at some criti- 
cal value of Uc and remains zero for all values oi U > Uc- 
The probability amplitude B does not vanish above Uc- 
Owing to this contribution, we have that for U > Uc 
the hopping of electrons from site j to a nearest neigh- 
bor is not forbidden, although restricted by the fact that 
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A = 0. This is consistent with the fact that the double 
occupancy remains finite for U > Uc and vanishes only 
in the limit U — > ooS2,. 

The phase diagram in the plane T-V is shown in 
Figs.[T21(top) and US (bottom). For U = 15 (see Fig. US 
(top)) there is no metallic phase and we have a criti- 
cal temperature Tc{V) where a transition from insulat- 
ing to charge ordered state is observed. The transition is 
first order up to Tc — 0.95, then becomes continuous. A 
re-entrant temperature is observed with the same char- 
acteristics previously discussed. For U = 8 (see Fig. El 
(bottom)) we have two critical temperatures, Td and Tc2, 
which characterize the MIT transition and the insulator- 
charge order transition, respectively. Also in this case a 
re-entrant temperature is observed in the latter transi- 
tion. 



V. DOUBLE OCCUPANCY, KINETIC AND 
INTERNAL ENERGIES 

In contrast to the local interaction U, a positive in- 
tersite interaction V favors the double occupancy D = 



W Ei(^*T(i)^*i(0) = f ~ ^22- This happens through the 
mechanism which favors the formation of the checker- 
board type charge ordered state as discussed above. The 
cost in energy U of having a double occupied site is partly 
balanced by the cost in energy V related to the pres- 
ence of two particles on nearest-neighbor sites. On the 
contrary, an attractive intersite Coulomb interaction V 
leads to a decrement of the double occupancy as the one- 
particle-per-site type of charge ordering favors the pres- 
ence of single occupied sites with respect to double occu- 
pied ones. 

A. One-dimensional system 

In Fig. El the double occupancy D is shown as a func- 
tion of the particle density for U ^ 2 and several values 
of the intersite Coulomb potential up to the critical value 
Vc [for U ~ 2 we have Vc « 2.47]. The observed features 
agree with the expectation that the double occupancy 
increases with the intersite potential. At half filling the 
double occupancy exhibits a discontinuity at the critical 
value Vc [for [/ = 4 we have Vc ~ 4.978] (see Fig. [T^ . 
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FIG. 14: The double occupancy versus the particle density at 
T = 0.01, (7 = 2 and several values of V . 
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FIG. 16: The kinetic energy K and the internal energy E 
are plotted versus the intersite Coulomb potential, for n = 1, 
T = 0.01 and various values of the onsite Coulomb potential. 



FIG. 15: The double occupancy at half filling versus V/Vc at 
T = 0.01 and several values of the onsite Coulomb potential. 

Again, this is a signal of a phase transition from the Mott 
insulator to a CO state. It is interesting to observe that 
at the middle point of the jump, the double occupancy 
takes the value 0.25 (the noninteracting one) for all val- 
ues of U . This result can also be inferred from Fig. 1 in 
Ref. Il^- It seems as the effect of V completely neutral- 
izes the effect of U and the system, at least for some local 
quantities, behaves as the noninteracting one. 

From the Hamiltonian H2.1|l we obtain the internal en- 
ergy per site Eh 

Eh ^^{H)= Adt{c'^{i)c\i)) + UD + dV{n°'(i)n{i)) 

(5.1) 

By increasing the intersite Coulomb potential the kinetic 
energy K = Adt{c"' {i)c^ {i)) decreases and the internal 
energy Eh increases, respectively, owing to the fact that 
the double occupancy increases. This is shown in Fig. 1161 



We can clearly see that the kinetic energy, as a function of 
the intersite Coulomb potential, shows a minimum at Vc 
when the transition is of the second order and develops a 
cusp when the transition is of the first order. The internal 
energy, instead, is almost insensible to the transition in 
the first case, but also develops a small cusp when the 
transition becomes of the first order (not shown) . 



B. Two-dimensional system 

In Fig. UZI the double occupancy is plotted versus the 
filling for different values of the intersite Coulomb poten- 
tial. The behavior is very similar to what has been found 
in the ID case. On increasing the intersite Coulomb po- 
tential the double occupancy behaves as that of a non- 
interacting system (i.e., tends to n^/4). 

The double occupancy as a function of the onsite 
Coulomb potential is studied in Figs. ^| (top) and El 
(bottom). In Fig. E| (top) we observe, by decreasing the 
onsite Coulomb potential, the transition from the insulat- 
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FIG. 17; The double occupancy versus the particle density at 
T = 0.01, (7 = 4 and several values of V . 
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FIG. 18: The double occupancy versus the onsite potential at 
T = 0.01, half filling and several values of V . 
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FIG. 19: The derivative of the double occupancy with respect 
to the onsite Coulomb potential is plotted versus the onsite 
potential at T = 0.01, half filling and several values of V . 



where V > & is studied in Fig. ^| (bottom); in accor- 
dance with the phase diagram shown in Fig. there is 
no metallic phase; by decreasing the onsite Coulomb po- 
tential, the double occupancy regularly increases until 
the charge ordered phase is reached. 

The phase transitions are marked by a discontinuity 
in the double occupancy. This is shown in Fig. 1191 where 
{dD/dU)n=i is plotted versus the onsite Coulomb poten- 
tial for different values of the intersite Coulomb potential. 
For V <Q (see Fig.^|(top)) we have two discontinuities, 
related to the MIT and insulator to charge ordered phase 
transitions. For V >Q (see Fig. ^] (bottom)) there is no 
metallic state and we observe only one discontinuity, re- 
lated to the latter type of transition. 



ing phase to the metallic phase and to the charge ordered 
phase. In the insulating state the double occupancy is 
quite depressed; at the MIT there is a discontinuity and 
the double occupancy rapidly increases by decreasing the 
onsite Coulomb potential until encounters a second dis- 
continuity when a CO phase comes into play. The case 



Accordingly to the behavior of the double occupancy, 
for positive values of the intersite Coulomb potential, the 
mobility of the electrons increases, signalling an overall 
tendency towards a charge density instability driven by 
the V term (see Fig. I2()|l . in agreement with the mean- 
field resulliS. 
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FIG. 20: The kinetic energy as a function of the particle den- 
sity for U — 4, T — 0.01 and various values of V . 



FIG. 22: The density of states at U 
V = 1. 
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FIG. 21: The Fermi surface for various values of V at f/ = 4 
and n = 0.7, 0.9. 



VI. THE FERMI SURFACE AND THE DOS 

With respect to what we have found in the sim- 
ple Hubbard modelM, the overall shape and bending 
of the Fermi surface of the system does not change on 
varying the intersite Coulomb potential, but, almost 
rigidly, its volume decreases on increasing the intersite 
Coulomb potential [see Figs.|2]- This can be explained 
as an isotropic increment of the available states in k- 
space and can be useful to describe quantitatively rather 
than qualitatively-"* (as the usual Hubbard model does) 
the ARPES determinations of the Fermi surface of the 
cuprate superconductors^^. 

As we can see from Figs. ED (top) and ED (bottom), 
for given values of the onsite Coulomb potential and the 
intersite Coulomb potential, there is a critical value of 
the doping ric where the Fermi surface is nested. At 
this critical value there is a crossing of the van Hove 
singularity and the Fermi level in the density of states 
(see Fig. 123 ■ 

This feature may be related, in the framework of 
the van Hove scenario, at the experimental results for 
the static susceptibilitySS*SI and specific heat^^i^^i'^°i'^^ 
in cuprate superconductors which exhibit well-definite 
peaks as functions of the filling. We have studied the 
value of the filling Uc, the result is shown in Fig. 1231 
where ric is plotted as a function of the intersite Coulomb 
potential for various values of the onsite Coulomb poten- 
tial. As we can see, ric increases on increasing the inter- 
site Coulomb potential until a certain value of C/ w 8 is 
reached. Above this critical value of the onsite Coulomb 
potential the influence of the intersite Coulomb potential 
is almost zero. This is a clear indication that if we would 
like to explain the cuprate superconductors and their 
anomalous features by means of Hubbard-like models, 
we need to exploit the intermediate regime for the onsite 
coupling and the weak regime for the intersite one. Only 
in this region of the model parameters, ric is in qualitative 
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FIG. 23: The critical value ric of the doping is shown as a 
function of V for various values of U. 



FIG. 25: The specific heat versus the temperature at n = 0.7, 
V — 1.5 and several values of U. 




FIG. 24: The specific heat versus the temperature at [/ = 4, 
n = 0.7 and several values of V. 



B. Two-dimensional case 

As regards the 2D case, it is worth noting that there is 
a crossing point when the curves are plotted for different 
values of the onsite Coulomb potential at fixed V (see 
Fig l25|l . but there is no crossing point when the onsite 
Coulomb potential is fixed and different values of the in- 
tersite Coulomb potential are considered. This can be 
seen as an indication of the fact that the two interaction 
terms act to different orders in a perturbation expansion 
for the double occupancy and the kinetic energy. In par- 
ticular, there is no value of temperature for which the 
first derivative of the double occupancy and the kinetic 
energy with respect to the intersite Coulomb potential 
does not depend on the intersite Coulomb potential. 



agreement with the experimental data&SLSiSSiZSiii. 



VIII. CONCLUSIONS 



VII. THE SPECIFIC HEAT 

The specific heat is shown in Fig. [21 for the one- 
dimensional case and in Fig l25l for the two-dimensional 
one. 



A. One-dimensional case 

In ID, the net effect of the intersite Coulomb repulsion 
is to reduce the splitting between the charge and spin 
energy scales with respect to what observed for the simple 
Hubbard model, with a resulting more pronounced single 
peak at intermediate temperatures (see Fig|] 



The extended Hubbard model, in one and two dimen- 
sions, has been studied by means of the Composite Oper- 
ator Method within the 2-pole approximation. According 
to the sign of the intersite Coulomb potential, transitions 
between the paramagnetic (metal and insulator) phase 
towards two kinds of charge ordered states (one-particle- 
per-site and checkerboard types) have been found. The 
rank of these transitions has been studied through the 
analysis of double occupancy, nearest-neighbor density- 
density correlator, kinetic and internal energies behavior. 
The evolution of the Fermi surface and of the nesting fill- 
ing has been tracked on changing the value of the intersite 
potential evidencing the range of parameters suitable to 
describe the experimental features of the cuprate super- 
conductors. The specific heat features on varying the 
intensity of the intersite potential has been contrasted to 
that found for the simple Hubbard model. 
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APPENDIX A: SPECTRAL DENSITY 
FUNCTIONS AND ENERGY SPECTRA 

Within the 2-pole approximation, the energy spectra 
Ei (k) and the spectral density functions cr^*^ (k) , appear- 
ing in Eq. 12.51 are given by 



^i(k) =i?(k) + Q(k) 

^2(k) = i?(k)-g(k) 



'12 



(k) 



111 
2 

hi 

2 



1 - 



g(k) 

2g(k) 

g(k) 

2g(k) 



"T-12(k) 

2g(k) 



(Al) 



(A2) 



III = 1 -~n/2 and /22 = ^^72 are the only non-zero entries 
of the normalization matrix; the parameters appearing in 
the previous equations are defined as 



A = Cli ~ C^^ 



(A7) 
(A8) 



We can see that the Green's function depends on the 
following set of parameters: ^ , Cii, (7^2, C^2 ' Pi X" ■ 
They can be computed as functions of the model param- 
eters and temperature and filling by the following set of 
coupled self-consistent equations 



"^l2(k) 
' 2g(k) 



I22 


\\ ^^^^ 1 

2g(k) 




I22 


[ 2g(k)J 





(A3) 



4^(k) 
4f(k) 

with the notation 
R = Ro + Ria{k) 

^ Jll-'22 ^ 

i?l - t4-{-*[P + ^22(1 - ")] + ^^[/llC^2"2 + ^22Ci"i]} 
-'ll-'22 

(A4) 

5(k) =90+ .gia(k) 

50 = -t/ + 7^ - «)A + - xc )] 

-'11J22 ^ 

2d 

91 = ^^Wl - - ^22) + V{l22C^l - IllC^2)] 



I11I2 

mi2(k) = mo + mia(k) 
mo = 2(it A 

mi = 2d[t{p-l22) + VC'l 



(A5) 



(A6) 



n = 2[1 - Cii - C22] 



(A9) 



^12 



n-2D 



(AlO) 



n - 2i? ^22 (^12 



where C^^ = (V'^(i)V'J («)) and D = ^ Ei('^T(i)"i(i)) = 
§ — C22 is the double occupancy. 

The correlation function C{i,j) = '^^'ii be 

computed in terms of the retarded propagator, or better 
in terms of the energy spectra and of the spectral density 
functions, by means of the following expressionSS 



n 



, AJh 



a(")(k)[l + r„(k)] (All) 



where and Hb are the volume of the unit cells in 
the direct and inverse spaces, respectively, and T„(k) — 
tanh {En{k)/2kBT). 
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